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1. Introduction 

1.1. Spectral Pollution in the Galerkin method. Let A be a self- 
adjoint operator acting on an infinite dimensional Hilbert space 7i, 
with a dense domain Dom{A). The spectrum of A, Spec{A), may 
be expressed as the union of the discrete spectrum consisting of all 
isolated eigenvalues of finite multiplicity, SpeC(jjsj.(y4), and the essential 
spectrum, where SpeCess(y4) := Spec(yl)\SpeC(jjgj.(74). In most standard 
situations the essential spectrum can be found analytically, but points 
in SpeC(jisc.(/l) are usually estimated by numerical procedures. 

The estimation of SpeC(jigj,(A) is often performed through subspaces 
C C Dom(yl) and corresponding truncations of A. Standard numerical 
techniques, such as the finite element method, aim at solving Galerkin 
approximate problems posed in weak form: 

, . find 7^ mG£ and AgM such that 

{Au,v) = X{u,v) \/v^C 

where C is finite dimensional. 

Backed by the Rayleigh-Ritz variational principle, when applicable, 
the Galerkin method represents a powerful tool in the analysis of spec- 
tral properties of linear operators. However, the Galerkin method is 
not foolproof, in general, the solutions of (P) might fail to provide re- 
liable information about the location of Spec (A) (see [3], [1], [7], [8], 

my 

The drawbacks in the Galerkin method are due in part to the so 
called spectral pollution phenomenon which we now describe. Let C 
Dom(/l) be a sequence of subspaces approaching ?i, as ra ^ oo (e.g. 
satisfying ([3]) below with p = 0, 1 only). Suppose we found 7^ n„ G £„ 
and A„ G M solutions of (P) with C = Cn, satisfying Xn —>■ fi and 
ll^nll^^^^n ^ w E Dom(y4) in the weak topology. By the approximating 
property of £„, we may obtain 

{Aw — fiw, v) = Vf G Dom{A), 

which appears to suggest that /i is in Spec(y4). Unfortunately, the 
latter conclusion is not ensured in general. Without further information 
about the structure of A (e.g. compactness properties), w might be 
G Dom{A), so convergent solutions of the approximate problem might 
produce "polluted" sequences A„ — > /i ^ Spec(74). 

The emergence of spurious eigenvalues in gaps of SpeC(,gg(/l) repre- 
sents a serious difficulty in applications such as elasticity theory and 
solid state physics (see [3] and [9]), as there is no universal recipe to 
detect or prevent them for a given operator A and sequence of approx- 
imate subspaces 
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1.2. Pollution-free strategies and quadratic methods. Spectral 
pollution is a consequence of the fact that in (P) we are truncating si- 
multaneously both u and v. Indeed, let 11 be the orthogonal projection 
onto C and 

^ . mx-A)v\\ 

-rrla;) := mm . 

^ ' O^vGC \\v\\ 

Then A G M satisfies (P) if, and only if, Fc{X) = 0. That is to say, 
there exists u E C such that (A — A)u ± C. As ||(A — is 
not guaranteed to be small, we have no indication whether A is close 
to Spec(y4) or not. 

This argument suggests that the correct quantity to look at is 

- A)v\\ 



Frix) := min 



As 



Fc{x)> inf M^^IM = = dist[x,Spec(A)], 

«eDom(A) ||n|| 

Fc{x) can be close to only when x is close to a point in the spectrum 
of A. 

In [8] , Davies and Plum considered a pollution-free strategy for find- 
ing Spec (A) based on computing the profile of Fc{x) for x G M. If 
£ C Dom(A2), 

FAxf = min ^ ^ „ ' ' ' = \mix - Af \ C]-^\\~^ 

^ ' . mix - A)'v\\ 

= mm — =: Grix). 

Oytv£C \\v\\ 

Therefore estimating Fc{x) reduces to computing eigenvalues of self- 
adjoint matrices depending on the parameter x G M. 

The approach developed in [8] relies heavily on being able to find 
accurately a matrix representation for n(x — A)^ f £ in terms of an 
orthonormal basis of C This is a drawback, for instance, if C is given 
by the finite element method, where an orthonormalisation of the basis 
will be numerically expensive. 

An alternative pollution-free method which is independent of the 
matrix representation of n(x — A)"^ \ L is also available and it may 
be obtained by considering the zeros of the function Gci^z) for 2 G C. 
Typically Guiz) and F\{z) only coincide at 2; G M. The (2dim£) 
zeros of the polynomial det (n(z — A)'^ \ are the zeros of Gc.{z) 
and, on the other hand, F[i{z) 7^ unless z is an eigenvalue of A, 
with corresponding eigenvector m G £, a very unlikely situation. It 
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is remarkable, however, that the non-real zeros of Gc{z) also provide 
reliable information about the location of Spec (A). 

This alternative procedure has been recently discussed in [7], [1] and 
[2], and it can be traced back to [1] and [10]. A central role is played 
by the problem 

. „ X find C G C such that 3n G £ with 

{Au.Av) ~2C{Au,v) + C{u,v) = Mv e C. 

It is readily seen that Gc{Cj = if, and only if, ^ is a solution of (Q). 
The philosophy of the method is to regard (Q), in place of (P), as an 
approximate spectral problem for operator A. 

The following universal non-pollution result justifies favouring (Q) 
over (P) (see [3, Theorem 2.6] or Theorem [3] below) : if C is a solution 
of (Q), then 

(2) dist[ReC,Spec(A)] < |ImC|. 

That is to say, ( can be close to M, only when it is also close to the 
spectrum of A. 

Problem (Q) gives rise to a matrix spectral problem quadratic in the 
spectral parameter. This added complication balances out with the 
reliability of the method expressed in the above result. 

Now, will a solution of (Q) ever be close to M? As for the Galerkin 
method, in general, additional conditions on a sequence of subspaces £„ 
are required for convergence. A precise statement reads as follows, see 
[2] or Theorem O below. Let A G Specjisj.(A) and n„ be the orthogonal 
projection onto Cn C Dom(A^). If 

WUnAm^u - XPu\\ Vp = 0,l,2, 

^ ' Vm G Dom(A) : Au = \u, 

then there exists Cn ^ C satisfying (Q) with C = Cn and £ = £„, such 
that |C„ — A| ^0. The above hypothesis is fulfilled immediately, for 
instance, if A is bounded and UnV v for all v ETi. 

The combination of these two results appears to provide a general 
pollution-free procedure for finding discrete eigenvalues of self-adjoint 
operators. Although this might seem too optimistic at the present 
moment, one of the advantages of this method lies in the fact that 
it is applicable without any special restriction upon the structure of 
Spec(/1). Moreover, the requirements on £„ are analogous to those 
needed in the Galerkin method. 

1.3. Stability of Quadratic Projection Methods. On the down- 
side, here we are confronted with a more difficult problem to solve. In 
general, the finite-dimensional eigenvalue problem associated to (Q) is 
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non-Hermitian. Accuracy, as well as stability of the method becomes 
a delicate matter. The main goal of the present note is to discuss how 
non-pollution and convergence of the method are affected, when the 
coefficients of problem (Q) are known only approximately. 

In Section [2] we will show that the non-pollution property remains 
stable in a sense which will be specified below. In Section |3] we will 
discuss stability of approximation. Note that a consistent formulation 
of (Q) only requires C Dom{A), see Remark m Under a suitable 
hypothesis on the subspaces £„, our Theorem [S] extends the analogous 
result of [2] by allowing n [Dom(yl) \ Dom(/l2)] ^ 0. In the final 
section we report on various numerical experiments performed on a 
simple example. 

2. Pollution-free Stability 

We devote this section to showing that, given error bounds in the 
computation of the coefficients of problem (Q), it is possible to con- 
trol errors in the pollution-free estimation of Spec (A) by the quadratic 
method described in Section [L2l 

Let us begin by fixing some notation. Below £„ is an n-dimensional 
subspace of Dom(y4) with basis {ei, . . . , e„}. This basis will always be 
normalised, ||ej|| = 1 for all j = 1, . . . ,n. When sufficiently clear from 
the context, we will suppress the sub-index and write £ = £„. 

For any u G C, u = ulci + ■ ■ ■ + ti^e„, from which we define the 
following norm on £, 

||m||o := H h \Un\^Y. 

Since £ is a finite dimensional space, there exists /9 > 0, such that 

(4) \\u\\ = {u,u)^ > [3\\u\\q \/ueC. 

If {ei,...,e„} is an orthonormal basis, then || ■ || = || • ||o- However 
when the basis is far from being orthonormal, /5 will be small. We will 
occasionally write u = {ui, . . . , n„) G C". 

Let matrices Aq, Ai and A2 in C*^^*^ be given entrywise by 

(5) [Aoljk = {Aej,Aek), [Ai]jk = {Aej,ek), [A2]jk = {ej,ek). 
Define the matrix polynomial M{z) G C"^" as 

(6) M{z) := Ao - 2zAi + 2 G C. 

Then C G C is a solution of (Q) if, and only if, det[M(C)] = 0. 

The stability results we establish below give a positive answer to the 
following question. Suppose we are only able to estimate the matrices 
Ap by Ap and the norm of the error \\Ap — Ap\\ < Bp, p = 0, 1, 2. Can 
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we recover information about the spectrum of A from the approximate 
problem 

(Q) find C e C : det [io - 2(^1 + (^^2] = 
with accuracy possibly depending upon e 



? 



The following preliminary result will be needed. 
Lemma 1. For 2 G C and 6 > 0, let 

J = [Re z — |Im z\ — S,Re z + |Im z\ + S], 
n = {{x- zY : X eR\J}. 

Then, 

(7) inf Re t; = 26\lm z\ + 5^. 

Proof. Let z = a + ib, a,b E M., then for any v E Q we have 

V = {x — a)^ — — 2b{x — a)i 
for some x G M\J . It is clear that Re f > and moreover 

inf Re t; = inf ix — a)^ — b'^ 

= 26\b\+6^ 

verifying ([7]). □ 

Theorem 2. Let A be a self-adjoint operator acting on a Hilbert space 
Ti, and C be an n- dimensional subspace o/Dom(y4). Let B be a singular 
n X n matrix. For any 2 G C, let M{z) and (3 be as in ([6]) and (j4]). 
Let G M with > \\M{z) - i?||c". If S > is such that (5^ + 
25|Im z\)P^ > a^, then 

(8) Spec(v4) n [Re z - |Im z| - 5, Re z + |Im z| + 5] 7^ 0. 

Proof. Let 5 > be as in the hypothesis and suppose the intersection 
([HD is empty. Using the spectral theorem and ([7]), we have for all m G £ 



Re {u^M{z)u) = Re ^ {{A - z)ej, {A - z)ek)ukUj 



jk=i 



= Re ^ {{A — z)ujej, {A — z)ukek) = Re {{A — z)u, {A — z)u) 
jk=i 

= [ Re{X- zf d{Exu,u) > {26\lmz\+6^)\\uf 
Jr 

> (25|Im z\ + 6^)(3^\\u\\l = {26\lm z\ + 6^)(3^\\u\\ln, 



STABILITY OF QUADRATIC PROJECTION METHODS 7 

where Ex is the spectral measure associated to A. It then follows from 
the Schwarz inequality that for any u G 

||M(^)u||c™ > {26\lmz\+d^)f3^\\u\\cn, 

so that the operator M{z) : C" — C" is invertible and 

\\M{z)-'\\cn < {{26\lmz\+6^)f3Y' <a-\ 

In particular ||M(2)~^ > ||M(2;) — -B||c", from which it follows that 
B is not singular. The result follows from the obtained contradiction. 

□ 

The next theorem is the main result of this section and it is an 
improvement on [3, Theorem 2.6]. 

Theorem 3. Let A be a self- adjoint operator acting on a Hilbert space 
H, and C be an n- dimensional subspace o/Dom(A). Define Aq, Ai and 
A2 as in Let Ap be n x n matrices, such that for Ep > 

\\Ap - ApWcn < Ep, p = 0,1,2. 

// the matrix Aq — 2(Ai + C^^2 is singular for some ( G C, then 

(9) Spec(A) n [Re C - 5, Re C + 5] ^ 
for 

5= V|ImCP + r^(ICP^2 + 2|Cki + eo). 

Proof With the notation of Theorem [21 take B = Aq - 2(^1 + C^^2- 
Since 

||M(C)-i?||c" < (|Cr^2 + 2|Cki + £o), 
(Q follows from Theorem [2l □ 

In particular, under the hypothesis above, 

(10) dist[Re C, Spec(yl)] < |Im C| + 

with e = A/|Cp£^2 + 2|Cki + £^0- If the basis {ei, . . . , e„} is orthonormal, 
then ([9]) and ffTOl) hold with (3 = 1. Note that the case eo = ^1 = ^2 = 
corresponds to [3, Theorem 2.6], see ([2]). 

3. Stability of Convergence in the Quadratic Method 

A consistent formulation of (Q) only requires C C Dom(A). How- 
ever, the available approximation results for the quadratic method (cf. 
[1] and [2]) impose the hypothesis C C Dom(/l^). In this section we 
show that, if C Dom(y4) approach reasonably well the eigenspace 
associated to an eigenvalue A G Spec^js^(yl), then solutions of (Q) will 
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converge to A in the large n limit, and the process remains stable under 
perturbation of the matrix coefficients of the polynomial M{z). 

Remark 4. Allowing the possibility of test spaces C ^ Dom{A'^) is 
only relevant when A is unbounded. If A is a differential operator of 
order 2m and the trial spaces are constructed using the finite element 
method, C C Dom(y4^) requires C^"^~^ conforming elements, while C C 
Dom(A) only requires C^^~^ conforming elements. The performance 
of the interpolation algorithm in the finite element method is usually 
compromised as m increases. 

Below we highlight explicitly the dependency on n of approximate 
subspaces and operators, so we denote matrices M{z) and Ap, corre- 
sponding to £„,, by M'^")(2r) and A^p \ respectively. We also assume 
throughout this section that the basis {ei, . . . , e^} of is orthonor- 
mal. In general we do not assume that C Cm whenever n < m. 
Strictly speaking we should denote the basis functions of Cn by {e^"''}. 
However we suppress this notation as no confusion shall arise. 

For u G Dom(74), the projection of u onto Cn is then given by 

n n 
k=l k=l 

Since {ei, . . . , e„} is orthonormal, we can isometrically identify £„, with 
C". 

Our key result assumes the following hypothesis on the sequence £„: 

, , Vp, g = 0, 1, and Vm G Dom(y4) : Au = Xu, 

^ ' II Y.%i{A^^nU, A'Jej)ej - ^0, as n ^ oo. 

Whenever £„ C Dom(yl^), (H) reduces to ([3]). Furthermore, if A is 
bounded and n„ converges strongly to the identity, then (H) holds 
true for all A G SpeC(jjgj,(/l). 

The following result is an improvement upon |2l Theorem 2.2]. 

Theorem 5. Let A be a self-adjoint operator on a Hilbert space. Sup- 
pose that the sequence of approximate subspaces Cn C Dom(yl) satisfy 
(H). Let A G SpeC(jisj,(74) and let d := dist[A, Spec(y4) \ {A}]. Given 
< 5 < d/A, there always exist N., 60,61,62 > ensuring the following. 
If n > N and the matrices Ap G C^^^ satisfy 

\\Ap — Ap II < 6p, p = 0,1, 2, 

then 

(a) we can always find ( E C with det[Ao — 2(Ai + C^^2] = and 
|C-A|<5, 
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(b) the set {/i G C : det[Ao — 2/_Mi + fi'^A2] = 0} does not intersect 
the annulus {w G C : 5 <\w — \\ < d/A}. 

The proof of this result will be given at the end of this section. It will 
be a consequence of various technical lemmas, in particular, suitable 
extensions of [21 Lemmas 5.1 and 5.3] and various regularity properties 
of Gc{z). 

We begin with the rigorous definition of the right hand side of ([1]) 
in the case C C Dom(y4). For 2; G C, let 

. ||M(z)v||cn 

Gc{z) := mm — - — . 

o^veC" ||v||c" 

If £ C Dom(A^), then Gc{z) coincides with the right hand side of ([1]). 
Below we will write Gn{z) = Gc„{z). 

Clearly ^^(C) = if, and only if, det M(^) = 0; so the solutions of 
problem (Q) are completely characterised as the zeros of Gc{z). It is 
readily seen that: 

(11) Gc{z) = \\[M{z)]-^\\^^ = least singular value of M(z). 

In fact Gciz)~^ is a continuous subharmonic function in the region 
G C : det M(2;) 7^ 0}, with singularities at the zeros of det[M(2;)] 
(see e.g. |1] or [21 Lemma 4.1]). This property will play a central role 
below. 

The statement of Theorem [5] will be obtained as a consequence of the 
fact that Gc{z) is small if, and only if, for small enough Sp > 0, z = ( 
is a solution of an approximate problem (Q). The following notion, 
which has recently become standard, will simplify considerably most 
of our arguments. Let 

A^{eo, ei, 62) ■.= {zeC : Gc{z) -{eo + 2si\z\+ eakH < 0}. 

This set is called the structured pseudospectrum of the matrix polyno- 
mial M{z), see [5]. 

The proof of the following fundamental property of the pseudospec- 
trum is a direct consequence of ffTTl) and [3 Lemma 2.1]. It clearly 
suggests how to verify the validity of (a) and (b) of Theorem [51 

Lemma 6. The complex number ( G A^(eo7 £^1, £^2) if, o-nd only if, 
det[io - 2Cii + C^^2] = for some Ap G C"''" satisfying \\Ap - Ap\\ < 
Ep, p = 0,1,2. 

Furthermore, cf. [6l Theorem 2.3], 
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Lemma 7. Let Q be a connected component of A'~{eo, 81,62), such that 
det M(/i) = for some /i G fi. // \\Ap — Ap\\ < Sp, there always exist 
QeVt such that det[io - 2C^i + (^M] = 0. 

We now establish two key relations between the large n limit of Gn{z) 
and dist[2;, Spec(/1)]^ in a neighbourhood of the discrete spectrum of 
A. 

Lemma 8. Let A G SpeC(jjgj,A. If the sequence of approximate subspaces 
Cn C Dom(A) satisfy (H), then 

lim G„(A) = 0. 

n— »oo 

Proof. Let u G Dom(74)\{0} be such that Au = Xu. Consider the 
vector u = ((ci, u), . . . , {cn, u)). We have 



[MW(A)u], = ^[M(A)],,(e,,u) 



n 

= ^{Aci, A{u, ej)ej) - 2X{Aei, {u, ej)ej) + A^(ei, {u, ej)ej) 
i=i 

= {Aci, Alinu) - 2\{Aei, n„M) + \^{ei, H^m), 

so that 

n 

||M(")(A)u||2„ = ^ |(An„M, Ae,) - 2A(n„M, Ae,-) + A'(n„u,e,)|2 

n 

= II ^^(An^M, Aej)ej — 2A(n„M, ^ej)ej + A^(n„M, ej)ej|p. 

The right hand side converges to zero by virtue of (H). Also ||u||c" 
||M||7^0as?T,^oo. Now, fix e > 0. Then, for all n large enough, 

G..(A)<"^^':'',f^"""<^- 

||u||c" 

As Gn{z) is non- negative and e is arbitrary the lemma follows. □ 

In general it is possible to construct examples where lim^^oo Gn{z) = 
for certain 2; ^ M, [1]. However, this is not possible for z in the vicinity 
of discrete eigenvalues of A. 

Lemma 9. Let A G Spec^iscl^) (md let d > be as in Theorem\^ As- 
sume that the sequence of approximate subspaces £„ C Dom(y4) satisfy 
(H). For all < d < d/4, there exist a constant < s < 1 such that 

(12) liminfG„(^) > s^^ for all 6 <\z- \\<d/A. 
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Proof. If Cn C Dom(A^), the result has been estabhshed in 
[21 Lemma 5.3]. We treat the more general case by considering ap- 
proximate subspace C Dom(A^) with orthonormal bases sufficiently 
close to Cn in the sense specified by (i)-(iii) below. 

As a first step, we recall the following standard result. For any 
V G Dom(y4), there exists a sequence w„ G Dom(A^) such that Vn ^ v 
and Avn — ^ Av. That is to say, Dom(y4^) is a core (in the operator 
sense) for A. 

Let 

c„ = max{l, max {\{A^ej, A''ek)\}\ ■ 

j,k=l...n;p,g=0,l 

Then it is always possible to find a set {ei, . . . , e„,} C Dom(A^) such 
that 

(i) {ei, . . . , e„} is orthonormal, 

(ii) ||ej- -e^ll < c^^e"", 

(iii) \{APe„A'^ek) - (A^e,, ^"4)1 < e-^ 

for j,k = 1, . . . , n and p,q = 0, 1. We may find ej by applying the 
Gram-Schmidt orthogonalisation procedure to vectors of Dom(y4^) suf- 
ficiently close to the ej. 

Let Cn '■= Span{ei, . . . ,e„} C Dom(y4^). In this proof, the symbol 
~ on top of matrices and operators denotes that they are constructed 
using C = Cn- Note that (iii) ensures the existence of complex numbers 
w^l such that \w^l\ < 1 and 

(A^e.^A'^ek) = {A^e,, A'^e,) + w^le'^ 

Property (iii) yields 

|[mW(z) - M^-H^)U = |2z[iS") - A^r^],, + [4") 
< 2\z\\{Ai,,ik) - {Aej,ek)\ + \{Aej,Aek) 
<{2\z\ + l)e~\ 

Thus, 

(13) ||M(")(2) - M(")(z)|| < {2\z\ + l)r2e-". 

Let u G Dom(y4) be such that Au = Xu. We next show that (ii) 
and the fact that (H) holds for Cn, ensures that (H) also holds for Cn- 



^0 \]k\ 



{Aej,Aeh 
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Indeed, 



k=l 



n 

< II Y,i^''^nU,A'^ek)ek - {Amr,u,A'^ek)ek\\ 

k=l 

n 

+ \\Y,{Amr,u,A^eu)ek-\P+M\ 

k=l 

n 

< II Y,{A''^nU,A''i,)ek - {Amr,u,A'^ek)ek\\ 

k=l 

n 

+ II J2(^''^nU,A''ek)ek - {Am^u,A'^ek)e4 

k=l 
n 

+ II J](APn„M, A''ek)ek - X^+'uW =T,+T2 + Tg. 



k=l 



Since satisfies condition (H), T3 — s> 0. We must show that the first 
two terms also converge to zero. Consider the first term, 



Ti = II 5];(n,g,)(A*'g„A'?4)(4-efc)|| 



jk=i 



< ||n|| ^ \{A%,A'iik)\\\ ek- ek\ 



jk=i 

n 

= \\u\\ J2 m''ej,A'^ek)+w^le-'')\\\ek-ek\\ 
jk=i 

n n 

< ||M||ne""^ 114 - Cfcll + ||m|| ^ \ {A^ej, A'^ek)\\\ek - ek\\. 

k=l jk=l 
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Using (ii) it is clear that Ti — > 0. For the second term we have, 

n 

T2 = \\Y, e,-) {A^e,, A''ek)ek - {u, e,) ^^6^)6^11 
jfc=i 

n 

= II - (]¥i~A^ej)efc|| 

3k=\ 

n 

< \\u\\ ^ IKA^Cj, A'^ek)ej - {A^Cj, A''ek)ej\\ 

jk=i 

n 

= \\u\\ J2 \\{{A^e„A'^ek) + w;^,'e-")e, - {A^e,, A'^e,)e,\\ 
jk=i 

n 

< ||M||n2e-" + ||m|| ^ KA^e,-, A%fc)|||g, - e,-||. 

jk=i 

Again, using (ii) it is clear that T2 0. This ensures that also 
satisfies (H) so (fT2l) is valid for Gn{z). 

The proof of ([12]) follows. Fix s > 0. Let v„ e C" such that ||v„|| = 1 
and 

||M(")(^)v„|| <G„(^)+5. 
Then, by virtue of (IT^ . 

Gn{z) < ||M("H^)Vn|| 

< ||[M("H^) - MW(^)]v„|| + ||mW(2)v„|| 

< (2|z| + l)ne-" + G„(z)+£. 

As this happens for all e > 0, the fact that Gn{z) satisfies f|T2|) implies 
that also Gn{z) satisfies this inequality. □ 

Proof of Theorem Let < s < 1 be as in Lemma M and d be as in 
the hypothesis of the theorem. By virtue of Lemmas M and [HI there 
exists > such that, Gn{\) < s6'^ and 

(14) Gn{z) > whenever 6 < \z - X\ < d/4, 

for all n > N. The subharmonicity of Gn{z)~^ ensures that the only 
local minima of Gn{z) are those points where the function vanishes. 
Thus, for all n > N, there always exists Cn £ C satisfying 



(15) 



ICn -\\<S and GniCn) = 0. 
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Let eo,ei,e2 > be small enough such that 

eo + 2ei\z\ + e2\z\'^ < s6'^ , \z-X\<d/4. 
Suppose that n > N. Then, by ffUl) . 

Gniz) - {eo + 2ei\z\ + esl^H > 

for all 5 <\z- \\ < d/A, so A^"{eo, Si, 62) fl {5 < |^ - A| < d/A} = 0. 
This, along with LemmaEl ensures (b). On the other hand, by virtue of 
(^Sii, A^"{eo,ei,e2)n{\z - X\ < 6} ^ 0. Thus Lemma [7] yield (a). □ 

4. Case Study 

Finite rank perturbations of multiplication operators have been con- 
sidered previously in connection with spectral pollution (see [8] , [7j and 
[2]) due to their simple structure. In this final section we report on var- 
ious numerical experiments we have performed on a model operator of 
this type. 

Let efc(x) = (27r)~2e*'''^ and 

..Jo for - TT < X < 0, 
"^^^ ~ \ 1 for < X < TT. 

In this section we assume that 

7^ = L2[-7r,7r], 

£2n+i = span{e_„(x), . . . ,e„(x)} and 

A(j){x) = a{x)(j){x) + {(j), eo)eo{x), (j) eH. 

Operator A is bounded and self-adjoint in Ti. Moreover, the spec- 
trum A is found explicitly. Since A is a rank one perturbation of the 
multiplication operator by the symbol a, Weyl's Theorem ensures that 
SpeCggg(y4) = Range(a) = {0, 1}. On the other hand, the isolated eigen- 
values of finite multiplicity of A are the solutions of 
((A — a)~^eo,eo) =!,[!]. A straightforward calculation reveals the 
two solutions = 1 ± which comprise the discrete spectrum of 

A. The eigenvalue A~ is inside the gap (0, 1) of the essential spectrum. 

As the symbol a(x) is discontinuous, the Fourier basis {cfc} is not a 
good choice for approximating A_ using the Galerkin method. Indeed, 
the solutions of (P) pollute the whole interval [0, 1] as the dimension 
of C2n+i increases. Let us test the quadratic method described in the 
preceding sections in this very simple model. 

Since A is bounded and 112^+10 — ^ 4> for all G 7i, condition (H) 
of Section [3] is satisfied. Thus, by virtue of Theorem [5l both discrete 
eigenvalues are approached by solutions of (Q) as n — >■ 00, free from 
spectral pollution according to Theorem [31 
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All the calculation described in this section were carried out using 
the computer package MATLAB. Fully functional m-codes are available 
at the web page [TT] . 

We compute the exact solutions of (Q), by finding the C G C such 
that det M^'^^^'^\() = 0. The matrix coefficients Ap^""*"^^ may be found 
explicitly using ([5]). They are sparse and Hermitian with entries either 
purely real or purely imaginary. The errors in solving (Q) are negligible 
for n of reasonable size (< 1000). 

In order to test the results established in the previous sections, we 
force large errors in the matrix entries, and compute the corresponding 
"perturbed" solution of the problem (Q). For simplicity, we fix Eq = 
ei = 62 = e. 

Let 



(16) [Ap],k = [Ap],k + 



2n + l 



-a 



(p) 

jk 



where a^} are random variables sampled from the unit disk {\z\ < 1} 
with additional constraints specified below. Then 



2n+l 



2n+l 

E£ (p) 
k=i 



2 2n+l 2n+l 

^ ' k=l jk=l 

2 2n+l 2n+l 2n+l 

£ E 1-^1' E 1 £ E I"'! 



2 



fc=l jk=l k=l 

SO II — Apll < e. Moreover this bound is sharp. Indeed, the matrix 
T such that [T]jk = for all 1 < j, < 2n + 1, satisfies 

||rf = ||T*r|| = £:||r||. 

We consider two types of restrictions on the random variable aj^^ 
On the one hand. Theorem [3] covers the general situation of moving all 
entries of Ap along randomly chosen directions in the complex plane. 

Thus, we perform unstructured perturbations by allowing all a^] 7^ 0. 
On the other hand, however, in order to reproduce the effect made by 
rounding errors in the estimation of the entries, we perform non-zero- 
Hermitian perturbations by imposing the condition: 



(p) _ 



0_ if [Ap\,k = 0, 
«SV0 if[Ap],,^0. 
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In Figure [H we depict the exact solutions of (Q) for n = 50. Ac- 
cording to ([2]), the points which are in the vicinity of the real axis are 
necessarily close to the spectrum. 

Figure [21 on the other hand, depicts the solutions of (Q) correspond- 
ing to 100 different random perturbations. Each of the graphs were 
constructed by prescribing a different constraint on the random vari- 
ables. Here n = 50 and e = 10"^. From the perturbed solutions one 
can identify Spec(74) less accurately but, once again, without pollution 
by virtue of Theorem [31 The correction 6 of Theorem [3l will depend on 
( and e, but notably not on n. Furthermore, Theorem [5] ensures that 
the clouds observed in Figure [2] will cluster near to each of the exact 
solutions of (Q) as e 0. 

Figures [3}[1| show the outcome of running Monte Carlo simulations 
in this model. We fix again e = 10~^. These pictures have been con- 
structed in the following manner. For each fixed n, we have found (~, 
the closest point to the eigenvalue A~ such that det M*^^"'+^)(^~) = 0. 
Then we have performed 20 constrained perturbations and averaged 
the solutions of (Q) which are closest to We know that solutions of 
the approximate problems close to (~ always exist, as a consequence 
of Theorem [5l Denote these averages by Cn'~ ^ind Q~ for the unstruc- 
tured and non-zero-Hermitian cases respectively. In Figure [3] we depict 
|Im (~\, |Im Cn'~l s-iid |Im Cn~\ n = 5 : 10 : 100. Similarly in 
Figure [H we depict |A~ — Re CrTl? |A~ — Re and |A~ — Re Cn~\- 

Figure [3] provides clear numerical evidence that the convergence of 
the quadratic method applied to this simple model is not lost even 
when the perturbations are large in modulus. Figures [H suggests that 
structured perturbations are considerably superior to the unstructured 
ones, in the test |A~ — Re C;7|- 

By combining Figure [3] and Theorem [3l we immediately predict a rate 
of convergence of |A~ — Re (~\ = o{n~^) for r f» 1/2. It is remarkable, 
however, that Figure [H strongly suggests an actual exponent of r ^ 1 
for this rate of convergence. An explanation of this phenomenon is 
linked to the fact that A~ is an isolated point of the spectrum, see [101 
Section 2]. We will be reporting on this issue elsewhere. 
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Figure 1. Exact solutions to (Q) for n = 50. 
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Figure 2. Top: solutions to (Q) for 100 unstructured 
random perturbations. Bottom: solutions to (Q) for 100 
non-zero-Hermitian random perturbations. Here e = 
10-1 
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Figure 3. Error predicted by Theorem [3] in the approx- 
imation of A~. Here we depict |Im (unperturbed), 
|Im Cn~\ |Im Cn~\ for n = 5 : 10 : 100. We average 
the two perturbed solutions of (Q) over a sample of 20 
problems with e = 10~^. The scaling is log-log and the 
horizontal axis shows 2n + 1. 




Figure 4. Actual error in the approximation of A~ us- 
ing the real part of solutions of (Q) and (Q) . Here we 
depict |A~ — Re (~\ (unperturbed), |A~ — Re Cn~\ 
|A~ — Re Cn~\ for ?^ = 5 : 10 : 100. We average the two 
perturbed solutions of (Q) over a sample of 20 problems 
with e = 10^^. The scaling is log-log and the horizontal 
axis shows 2n -|- 1. 



